//#  define CUSTOM_CUDA_SAFE_CALL(call) {                                      \
//    cudaError err = call;                                                    \
//    if( cudaSuccess != err) {                                                \
//        fprintf(stderr, "Cuda error in file '%s' in line %i : %s.\n",        \
//                __FILE__, __LINE__, cudaGetErrorString( err) );              \
//		getchar();                                                           \
//		exit(EXIT_FAILURE);                                                  \
//		}																	 \
//	}

/*
 * This file contains all methods related to image preprocessing as well as 
 * creating a 1d array from the image and segmenting it
*/

// Create a histogram of the image
//__global__ void createHistogramCUDA(unsigned char *d_imageArray, short *lut )
__global__ void createHistogramCUDA(unsigned char *d_imageArray, int *lut )
{
	int idx = (blockIdx.x * blockDim.x) + threadIdx.x;

	int moy = d_imageArray[idx];
	// We require an atomicAdd instruction because of the problem of race conditions.
	// This method of creating a histogram does not fully utilse the potential of the GPU as every thread is calling an atomic
	// add instruction which means the required value in global memory must be read, modified and rewritten back to global memory which
	// has a large latency.
	atomicAdd(&(lut[moy]), 1);
}


//__global__ void calEnhancementCUDA(unsigned char *d_imageArray, unsigned char *d_I1AndI2, float *d_M_C)
//__global__ void calLutCUDA(short *d_imageArray, unsigned char *d_I1AndI2, float *d_M_C)
__global__ void calculateI1AndI2(unsigned char *d_imageArray, unsigned char *pMinResults, unsigned char *pMaxResults)
{
	// Declare arrays to be in shared memory.
	// If THREADS_PER_BLOCK=32
	// THREADS_PER_BLOCK elements * (1 byte per element) * 2 = 64B
	// so each SM can have 16KB / (64B / 1024) = 256 blocks
	__shared__ unsigned char min[THREADS_PER_BLOCK];
	__shared__ unsigned char max[THREADS_PER_BLOCK];


	unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x*(blockDim.x*2) + threadIdx.x;

    unsigned char val1 = d_imageArray[i];
	unsigned char val2 = d_imageArray[i + THREADS_PER_BLOCK];
	min[tid] = (val1 < val2) ? val1 : val2;
	max[tid] = (val1 > val2) ? val1 : val2;

    __syncthreads();

    // do reduction in shared mem
    //for(unsigned int s=blockDim.x/2; s>0; s>>=1)
    for(unsigned int s=blockDim.x/2; s>32; s>>=1)
    {
        if (tid < s) 
        {
			unsigned char temp = min[tid + s];

			if (temp < min[tid])
				min[tid] = temp;

			temp = max[tid + s];

			if (temp > max[tid])
					max[tid] = temp;

        }
        __syncthreads();
    }

	// unroll the last warp
	if (tid < 32)
    {
		{		
			min[tid] = (min[tid] < min[tid + 32]) ? min[tid] : min[tid + 32];
			max[tid] = (max[tid] > max[tid + 32]) ? max[tid] : max[tid + 32];
		}
		{		
			min[tid] = (min[tid] < min[tid + 16]) ? min[tid] : min[tid + 16];
			max[tid] = (max[tid] > max[tid + 16]) ? max[tid] : max[tid + 16];
		}
		{		
			min[tid] = (min[tid] < min[tid + 8]) ? min[tid] : min[tid + 8];
			max[tid] = (max[tid] > max[tid + 8]) ? max[tid] : max[tid + 8];
		}
		{		
			min[tid] = (min[tid] < min[tid + 4]) ? min[tid] : min[tid + 4];
			max[tid] = (max[tid] > max[tid + 4]) ? max[tid] : max[tid + 4];
		}
		{		
			min[tid] = (min[tid] < min[tid + 2]) ? min[tid] : min[tid + 2];
			max[tid] = (max[tid] > max[tid + 2]) ? max[tid] : max[tid + 2];
		}
		{		
			min[tid] = (min[tid] < min[tid + 1]) ? min[tid] : min[tid + 1];
			max[tid] = (max[tid] > max[tid + 1]) ? max[tid] : max[tid + 1];
		}
    }

	// write result for this block to global mem
    if (tid == 0)
	{
		pMinResults[blockIdx.x] = min[0];
		pMaxResults[blockIdx.x] = max[0];
	}

	/* From a previous version of calculateI1AndI2 (calEnhancementCUDA)

		////// Set only one single thread to carry out the gradient and intercept calculation.
		////if (idx == 640*480-100000)
		////{
		//	// Gradient
		//	d_M_C[0] = 255 / ((float)d_I1AndI2[1] - (float)d_I1AndI2[0]); // 255 / (I2 - I1)
		//	// Intercept
		//	d_M_C[1] = -d_M_C[0] * (float)d_1AndI2[0];
		////}
	*/
}


__global__ void calculateI1AndI2Stage2(unsigned char *pMinResults, unsigned char *pMaxResults, unsigned char *d_smallestValues, unsigned char *d_largestValues)
{
	__shared__ unsigned char min[THREADS_PER_BLOCK];
	__shared__ unsigned char max[THREADS_PER_BLOCK];


	unsigned int tid = threadIdx.x;
    unsigned int i = (blockIdx.x*blockDim.x) + threadIdx.x;

	//unsigned char minVal1 = pMinResults[i];
	//unsigned char minVal2 = pMinResults[i + THREADS_PER_BLOCK];
	//min[tid] = (minVal1 < minVal2) ? minVal1 : minVal2;

	//unsigned char maxVal1 = pMaxResults[i];
	//unsigned char maxVal2 = pMaxResults[i + THREADS_PER_BLOCK];
	//max[tid] = (maxVal1 > maxVal2) ? maxVal1 : maxVal2;

	min[tid] = pMinResults[i];
	max[tid] = pMaxResults[i];

	__syncthreads();

	for(unsigned int s=blockDim.x/2; s>32; s>>=1)
	{
		if ((tid < s) && ((i+s) < (MAX_DATA_SIZE / THREADS_PER_BLOCK / 2)))
		{
			unsigned char temp = min[tid + s];

			if (temp < min[tid])
				min[tid] = temp;

			temp = max[tid + s];

			if (temp > max[tid])
					max[tid] = temp;

		}
		__syncthreads();
	}

	// unroll the last warp
	if (tid < 32)
	{
		{		
			min[tid] = (min[tid] < min[tid + 32]) ? min[tid] : min[tid + 32];
			max[tid] = (max[tid] > max[tid + 32]) ? max[tid] : max[tid + 32];
		}
		{		
			min[tid] = (min[tid] < min[tid + 16]) ? min[tid] : min[tid + 16];
			max[tid] = (max[tid] > max[tid + 16]) ? max[tid] : max[tid + 16];
		}
		{		
			min[tid] = (min[tid] < min[tid + 8]) ? min[tid] : min[tid + 8];
			max[tid] = (max[tid] > max[tid + 8]) ? max[tid] : max[tid + 8];
		}
		{		
			min[tid] = (min[tid] < min[tid + 4]) ? min[tid] : min[tid + 4];
			max[tid] = (max[tid] > max[tid + 4]) ? max[tid] : max[tid + 4];
		}
		{		
			min[tid] = (min[tid] < min[tid + 2]) ? min[tid] : min[tid + 2];
			max[tid] = (max[tid] > max[tid + 2]) ? max[tid] : max[tid + 2];
		}
		{		
			min[tid] = (min[tid] < min[tid + 1]) ? min[tid] : min[tid + 1];
			max[tid] = (max[tid] > max[tid + 1]) ? max[tid] : max[tid + 1];
		}
	}

	// write result for this block to global mem
	if (tid == 0)
	{
		d_smallestValues[blockIdx.x] = min[0];
		d_largestValues[blockIdx.x] = max[0];
	}
}

/* Doesnt dynamically allocate shared memory - only works for THREADS_PER_BLOCK = 256

	//__global__ void calculateI1AndI2Final(unsigned char *pMinResults, unsigned char *pMaxResults, float *d_M_C)
	//{
	//	extern __shared__ unsigned char s_min[];
	//	extern __shared__ unsigned char s_max[];
	//
	//	//printf("pMinResults[0]: %d\n", pMinResults[0]);
	//	//printf("pMinResults[1]: %d\n", pMinResults[1]);
	//	//printf("pMinResults[2]: %d\n", pMinResults[2]);
	//
	//	unsigned char I1=255;
	//	s_min[0] = pMinResults[0];
	//	s_min[1] = pMinResults[1];
	//	s_min[2] = pMinResults[2];
	//
	//	//printf("s_min[0]: %d\n", s_min[0]);
	//	//printf("s_min[1]: %d\n", s_min[1]);
	//	//printf("s_min[2]: %d\n", s_min[2]);
	//
	//	unsigned char I2=0;
	//	s_max[0] = pMaxResults[0];
	//	s_max[1] = pMaxResults[1];
	//	s_max[2] = pMaxResults[2];
	//
	//
	//	for(unsigned int i=0; i<3; i++)
	//	{
	//		if(s_min[i] < I1)
	//			I1 = s_min[i];
	//
	//		if(s_max[i] > I2)
	//			I2 = s_max[i];
	//	}
	//
	//	//printf("I1: %d", I1);
	//	//printf("I2: %d", I2);
	//
	//	// Gradient
	//	d_M_C[0] = 255 / ((float)I2 - (float)I1); // 255 / (I2 - I1)
	//	// Intercept
	//	d_M_C[1] = -d_M_C[0] * ((float)I1);
	//}
*/

__global__ void calculateI1AndI2Final(unsigned char *pMinResults, unsigned char *pMaxResults, unsigned char minBlocksRequiredAtStage2, float *d_M_C)
{
	extern __shared__ unsigned char s_minAndMax[];


	unsigned char* minArray = &s_minAndMax[0];
	unsigned char* maxArray = &s_minAndMax[minBlocksRequiredAtStage2];

	unsigned char I1=255;
	for(unsigned int i=0; i<minBlocksRequiredAtStage2; i++)
		minArray[i] = pMinResults[i];

	unsigned char I2=0;
	for(unsigned int i=0; i<minBlocksRequiredAtStage2; i++)
		maxArray[i] = pMaxResults[i];

	for(unsigned int i=0; i<minBlocksRequiredAtStage2; i++)
	{
		if(minArray[i] < I1)
			I1 = minArray[i];

		if(maxArray[i] > I2)
			I2 = maxArray[i];
	}

	// Gradient
	d_M_C[0] = 255 / ((float)I2 - (float)I1); // 255 / (I2 - I1)
	// Intercept
	d_M_C[1] = -d_M_C[0] * ((float)I1);
}

//__global__ void test(float I_PixelValue, float gradient, float intercept, float *outPixelValue)
//{
//
//
//	unsigned char pixelValue = (unsigned char)I_PixelValue;
//
//	//printf("pixelValue: %d\n", pixelValue);
//	//printf("gradient: %f\n", gradient);
//	//printf("intercept: %f\n", intercept);
//	//printf("new pixelValue: %f\n", (gradient * pixelValue));
//	//printf("new pixelValue: %f\n", ((gradient * pixelValue) + intercept));
//
//	////outPixelValue[0] = ((gradient * pixelValue) + intercept);
//	outPixelValue[0] = (unsigned char) ((gradient * pixelValue) + intercept);
//}

__global__ void enhanceContrastCUDA(unsigned char *d_imageArray, float *d_M_C)
{
	int idx = (blockIdx.x * blockDim.x) + threadIdx.x;

	unsigned char pixelValue = d_imageArray[idx];
	float gradient = d_M_C[0];
	float intercept = d_M_C[1];

	if (pixelValue < (-intercept)/gradient)
		d_imageArray[idx] = 0;
	else if (pixelValue > (255 - intercept)/gradient)
		d_imageArray[idx] = 255;
	else
	{
		//if(idx==235826 || idx==224957)
		//{
		//	//printf("idx: %d\n", idx);
		//	//printf("idx value: %d\n", d_imageArray[idx]);
		//	//printf("gradient: %f\n", gradient);
		//	//printf("intercept: %f\n", intercept);
		//	//printf("new pixelValue: %f\n", (gradient * pixelValue));
		//	//printf("new pixelValue: %f\n", ((gradient * pixelValue) + intercept));

		//	//pixelValue: 217
		//	//gradient: 1.268657
		//	//intercept: -20.298508
		//	//new pixelValue: 255
		//}

		//if(idx==235826 || idx==224957) {
		//	d_imageArray[idx] = 255;
		//} 
		//else
			d_imageArray[idx] = (unsigned char) ((gradient * (float)pixelValue) + intercept);
	}
}

/* Not utilised

	//__global__ void enhanceContrastCUDAUsingLut(unsigned char *d_imageArray, float *d_M_C, short* lut)
	//{
	//	int idx = (blockIdx.x * blockDim.x) + threadIdx.x;
	//
	//	int pixelValue = d_imageArray[idx];
	//
	//	if (pixelValue < -d_M_C[1]/d_M_C[0])
	//		d_imageArray[idx] = 0;
	//	else if (pixelValue > (255 - d_M_C[1])/d_M_C[0])
	//		d_imageArray[idx] = 255;
	//	else
	//		d_imageArray[idx] = (int) (lut[pixelValue]);
	//}
*/

__global__ void lowPassFilterCUDATexture(unsigned char *r_imageArray)
{

	int idx = (blockIdx.x * blockDim.x) + threadIdx.x;
	
		short sum = 0;
		short count = 0;

		if (idx >= 640 )							{ sum += tex1Dfetch(texRef, idx-640);		count++; }		// Pixel above
		if (idx <= 306559 )							{ sum += tex1Dfetch(texRef, idx+640);		count++; }		// Pixel below
		if (idx % 640 != 0 )						{ sum += tex1Dfetch(texRef, idx-1); 		count++; }		// Pixel left
		if (idx % 640 != 639 )						{ sum += tex1Dfetch(texRef, idx+1); 		count++; }		// Pixel right
		if ((idx >= 640) && (idx % 640 != 639))		{ sum += tex1Dfetch(texRef, idx-640+1); 	count++; }		// Top right
		if ((idx <= 306559) && (idx % 640 != 639))	{ sum += tex1Dfetch(texRef, idx+640+1); 	count++; }		// Bottom right
		if ((idx >= 640) && (idx % 640 != 0))		{ sum += tex1Dfetch(texRef, idx-640-1); 	count++; }		// Top left
		if ((idx <= 306559) && (idx % 640 != 0))	{ sum += tex1Dfetch(texRef, idx+640-1); 	count++; }		// Bottom left

		sum += tex1Dfetch(texRef, idx); count++;

		r_imageArray[idx] = (unsigned char) (sum/count);
}


__global__ void lowPassFilterCUDA(unsigned char *d_imageArray, unsigned char *r_imageArray)
{

	int idx = (blockIdx.x * blockDim.x) + threadIdx.x;
	
	//if ((idx >= 640) && (idx <= 306559) && (idx % 640 != 0) && (idx % 640 != 639))	// Ignore all the border pixels	
	//{
		// Mask size of 9.
		short sum = 0;
		short count = 0;

		// 640 = WIDTH
		// 306559 = MAX_DATA_SIZE - WIDTH - 1

		if (idx >= 640 )							{ sum += d_imageArray[idx-640];		count++; }		// Pixel above
		if (idx <= 306559 )							{ sum += d_imageArray[idx+640];		count++; }		// Pixel below
		if (idx % 640 != 0 )						{ sum += d_imageArray[idx-1];		count++; }		// Pixel left
		if (idx % 640 != 639 )						{ sum += d_imageArray[idx+1];		count++; }		// Pixel right
		if ((idx >= 640) && (idx % 640 != 639))		{ sum += d_imageArray[idx-640+1];	count++; }		// Top right
		if ((idx <= 306559) && (idx % 640 != 639))	{ sum += d_imageArray[idx+640+1];	count++; }		// Bottom right
		if ((idx >= 640) && (idx % 640 != 0))		{ sum += d_imageArray[idx-640-1];	count++; }		// Top left
		if ((idx <= 306559) && (idx % 640 != 0))	{ sum += d_imageArray[idx+640-1];	count++; }		// Bottom left

		sum += d_imageArray[idx]; count++;


		//if (idx == 642)
		//{
		//	printf("Top Left    : %d\n", (d_imageArray[idx-640-1]));
		//	printf("Top         : %d\n", (d_imageArray[idx-640]));
		//	printf("Top Right   : %d\n", (d_imageArray[idx-640+1]));
		//	printf("    Left    : %d\n", (d_imageArray[idx-1]));
		//	printf("    Middle  : %d\n", (d_imageArray[idx]));
		//	printf("    Right   : %d\n", (d_imageArray[idx+1]));
		//	printf("Bottom Left : %d\n", (d_imageArray[idx+640-1]));
		//	printf("Bottom      : %d\n", (d_imageArray[idx+640]));
		//	printf("Bottom Right: %d\n", (d_imageArray[idx+640+1]));
		//	printf("Sum: %d\n", sum);
		//	printf("Count: %d\n", count);
		//	printf("r_imageArray[idx]: %d\n", sum/count);
		//	printf("r_imageArray[idx]: %d\n", (unsigned char) (sum/count));
		//	printf("\n");

		//	r_imageArray[idx] = 0;
		//}
		//else
		//{
			r_imageArray[idx] = (unsigned char) (sum/count);
		//}
	//}
}

/*

Isodata algorithm - This iterative technique for choosing a threshold was developed by Ridler and Calvard . The histogram is initially segmented into two parts using a starting 
threshold value such as 0 = 2B-1, half the maximum dynamic range. The sample mean (mf,0) of the gray values associated with the foreground pixels and the sample mean (mb,0) of the 
gray values associated with the background pixels are computed. A new threshold value 1 is now computed as the average of these two sample means. The process is repeated, based 
upon the new threshold, until the threshold value does not change any more. In formula:

*/
__global__ void calculateThresholdCUDA(unsigned char *d_imageArray, unsigned char d_threshold, unsigned int *d_greaterThanSum , unsigned int *d_greaterThanCount, unsigned int *d_lessThanSum)
{
	// 3072KB of shared memory.
	__shared__ unsigned int greaterThanThreshold[THREADS_PER_BLOCK];
	__shared__ unsigned int greaterThanThresholdCount[THREADS_PER_BLOCK];
	__shared__ unsigned int lessThanThreshold[THREADS_PER_BLOCK];

	//__shared__ unsigned int greaterThanThreshold[32];
	//__shared__ unsigned int greaterThanThresholdCount[32];
	//__shared__ unsigned int lessThanThreshold[32];

	//printf("$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n");
	//printf("BLOCK ID: %d\n\n", blockIdx.x);

	unsigned char threshold = d_threshold;
	//printf("d_threshold: %d\n", threshold);

	unsigned int tid = threadIdx.x;
    unsigned int i = (blockIdx.x*blockDim.x) + threadIdx.x;

	if (d_imageArray[i] > threshold)
	{
		greaterThanThreshold[tid] = d_imageArray[i];
		//greaterThanThresholdCount[tid]++;
		greaterThanThresholdCount[tid] = 1;
        //printf("greaterThanThreshold[%d]: %d\n", tid, greaterThanThresholdCount[tid]);
		lessThanThreshold[tid] = 0;

		//printf("greaterThanThreshold[%d]: %d\n", tid, greaterThanThreshold[tid]);
		//printf("greaterThanThresholdCount[%d]: %d\n", tid, greaterThanThresholdCount[tid]);
		//printf("lessThanThreshold: %d\n", lessThanThreshold[0]);
	}
	else
	{
		lessThanThreshold[tid] = d_imageArray[i];

		greaterThanThreshold[tid] = 0;
		greaterThanThresholdCount[tid] = 0;
	}




	__syncthreads();

	

    for(unsigned int s=blockDim.x/2; s>0; s>>=1) 
    {
        if (tid < s) 
        {
			//if(tid == 0)
			//{
			//	printf("tid: %d\n", tid);
			//	printf("s: %d\n", s);
			//	printf("greaterThanThreshold[tid]: %d\n", greaterThanThreshold[tid]);
			//	printf("greaterThanThreshold[tid + s]: %d\n", greaterThanThreshold[tid + s]);
			//	printf("greaterThanThreshold[tid]: %d\n", greaterThanThresholdCount[tid]);
			//	printf("greaterThanThresholdCount[tid + s]: %d\n", greaterThanThresholdCount[tid + s]);
			//	printf("lessThanThreshold[tid]: %d\n", lessThanThreshold[tid]);
			//	printf("lessThanThreshold[tid + s]: %d\n", lessThanThreshold[tid + s]);
			//	printf("------------------\n");
			//}

			greaterThanThreshold[tid] += greaterThanThreshold[tid + s];
			greaterThanThresholdCount[tid] += greaterThanThresholdCount[tid + s];
			lessThanThreshold[tid] += lessThanThreshold[tid + s];
        }
        __syncthreads();
    }



    if (tid == 0)
	{
 		d_greaterThanSum[blockIdx.x] = greaterThanThreshold[0];
		d_greaterThanCount[blockIdx.x] = greaterThanThresholdCount[0];
		d_lessThanSum[blockIdx.x] = lessThanThreshold[0];
	}
}

__global__ void calculateThresholdCUDAStage2(unsigned int *d_greaterThanSum , unsigned int *d_greaterThanCount, unsigned int *d_lessThanSum, \
											 unsigned int *d_greaterThanSumFinalValues, unsigned int *d_greaterThanCountFinalValues, unsigned int *d_lessThanSumFinalValues)
{
	// 3072KB of shared memory.
	__shared__ unsigned int greaterThanThreshold[THREADS_PER_BLOCK];
	__shared__ unsigned int greaterThanThresholdCount[THREADS_PER_BLOCK];
	__shared__ unsigned int lessThanThreshold[THREADS_PER_BLOCK];

	unsigned int tid = threadIdx.x;
    unsigned int i = (blockIdx.x*blockDim.x) + threadIdx.x;

	greaterThanThreshold[tid] = d_greaterThanSum[i];
	greaterThanThresholdCount[tid] = d_greaterThanCount[i];
	lessThanThreshold[tid] = d_lessThanSum[i];


	__syncthreads();

    for(unsigned int s=blockDim.x/2; s>0; s>>=1) 
    {
        if ((tid < s) && ((i+s) < (MAX_DATA_SIZE / THREADS_PER_BLOCK)))
        {
			greaterThanThreshold[tid] += greaterThanThreshold[tid + s];
			greaterThanThresholdCount[tid] += greaterThanThresholdCount[tid + s];
			lessThanThreshold[tid] += lessThanThreshold[tid + s];
        }
        __syncthreads();
    }


    if (tid == 0)
	{
 		d_greaterThanSumFinalValues[blockIdx.x] = greaterThanThreshold[0];
		d_greaterThanCountFinalValues[blockIdx.x] = greaterThanThresholdCount[0];
		d_lessThanSumFinalValues[blockIdx.x] = lessThanThreshold[0];
	}
}

__global__ void calculateThresholdCUDAFinal(unsigned int *d_greaterThanSumFinalValues, unsigned int *d_greaterThanCountFinalValues, unsigned int *d_lessThanSumFinalValues, \
											unsigned int minBlocksRequiredForThresholdingAtStage2, unsigned char *newThresholdValue)
{
	extern __shared__ unsigned int s_greaterThan_Count_lessThan[];

	unsigned int* greaterThanSumArray = &s_greaterThan_Count_lessThan[0];
	unsigned int* greaterThanCountArray = &s_greaterThan_Count_lessThan[minBlocksRequiredForThresholdingAtStage2];
	unsigned int* lessThanArray = &s_greaterThan_Count_lessThan[minBlocksRequiredForThresholdingAtStage2 * 2];

	unsigned int backgroundSetValue = 0;
	for(unsigned int i=0; i<minBlocksRequiredForThresholdingAtStage2; i++)
		greaterThanSumArray[i] = d_greaterThanSumFinalValues[i];

	unsigned int backgroundSetCount = 0;
	for(unsigned int i=0; i<minBlocksRequiredForThresholdingAtStage2; i++)
		greaterThanCountArray[i] = d_greaterThanCountFinalValues[i];

	unsigned int objectSetValue = 0;
	for(unsigned int i=0; i<minBlocksRequiredForThresholdingAtStage2; i++)
		lessThanArray[i] = d_lessThanSumFinalValues[i];

	for(unsigned int i=0; i<minBlocksRequiredForThresholdingAtStage2; i++)
	{
		backgroundSetValue += greaterThanSumArray[i];
		backgroundSetCount += greaterThanCountArray[i];
		objectSetValue += lessThanArray[i];
	}

	//backgroundSetCount = (backgroundSetCount == 0) ? 1 : backgroundSetCount;
	unsigned int backgroundSetAverage = backgroundSetValue / backgroundSetCount;
	//printf("backgroundSetAverage: %d\n", backgroundSetAverage);
	unsigned int objectSetAverage = objectSetValue / (MAX_DATA_SIZE - backgroundSetCount);
	//printf("objectSetAverage: %d\n", objectSetAverage);

	*newThresholdValue = (unsigned char) ((objectSetAverage + backgroundSetAverage) / 2);
	//printf("newThresholdValue: %d\n",*newThresholdValue);
}

int iterativeCalculateThresholdCUDA(unsigned char *d_imageArray)
//int iterativeCalculateThresholdCUDA(unsigned char *imageArray)
{
	int thresholdFound;

	unsigned char thresholdValue, newThresholdValue;
	newThresholdValue = 128;  //Half the dynamic range
	thresholdFound = 0;

	int minMemoryRequiredForThresholding = MAX_DATA_SIZE / THREADS_PER_BLOCK;
	
	//unsigned int* h_greaterThanSum = (unsigned int*)malloc(sizeof(unsigned int) * minMemoryRequiredForThresholding);
	//unsigned int* h_greaterThanCount = (unsigned int*)malloc(sizeof(unsigned int) * minMemoryRequiredForThresholding);
	//unsigned int* h_lessThanSum = (unsigned int*)malloc(sizeof(unsigned int) * minMemoryRequiredForThresholding);

	unsigned int* d_greaterThanSum;
	unsigned int* d_greaterThanCount;
	unsigned int* d_lessThanSum;
    CUSTOM_CUDA_SAFE_CALL( cudaMalloc( (void **)&d_greaterThanSum, sizeof(unsigned int) * minMemoryRequiredForThresholding) );
    CUSTOM_CUDA_SAFE_CALL( cudaMalloc( (void **)&d_greaterThanCount, sizeof(unsigned int) * minMemoryRequiredForThresholding) );
	CUSTOM_CUDA_SAFE_CALL( cudaMalloc( (void **)&d_lessThanSum, sizeof(unsigned int) * minMemoryRequiredForThresholding) );


	int minBlocksRequiredForThresholdingAtStage2 = (int) ceil(((float)minMemoryRequiredForThresholding / THREADS_PER_BLOCK));
	

	//unsigned int* h_greaterThanSum = (unsigned int*)malloc(sizeof(unsigned int) * minBlocksRequiredForThresholdingAtStage2);
	//unsigned int* h_greaterThanCount = (unsigned int*)malloc(sizeof(unsigned int) * minBlocksRequiredForThresholdingAtStage2);
	//unsigned int* h_lessThanSum = (unsigned int*)malloc(sizeof(unsigned int) * minBlocksRequiredForThresholdingAtStage2);

	unsigned int* d_greaterThanSumS2;
	unsigned int* d_greaterThanCountS2;
	unsigned int* d_lessThanSumS2;
	CUSTOM_CUDA_SAFE_CALL( cudaMalloc( (void **)&d_greaterThanSumS2, sizeof(unsigned int) * minBlocksRequiredForThresholdingAtStage2) );
    CUSTOM_CUDA_SAFE_CALL( cudaMalloc( (void **)&d_greaterThanCountS2, sizeof(unsigned int) * minBlocksRequiredForThresholdingAtStage2) );
	CUSTOM_CUDA_SAFE_CALL( cudaMalloc( (void **)&d_lessThanSumS2, sizeof(unsigned int) * minBlocksRequiredForThresholdingAtStage2) );

	unsigned char* h_threshold = (unsigned char*) malloc(sizeof(unsigned char));
	unsigned char* d_threshold;
	CUSTOM_CUDA_SAFE_CALL( cudaMalloc( (void **)&d_threshold, sizeof(unsigned char)) );

	unsigned int sharedMemSize = 3 * minBlocksRequiredForThresholdingAtStage2 * sizeof(unsigned int);
	


	while(thresholdFound != 1)
	{
		//printf("\n****************\n");
		//printf("thresholdValue - Before: %d\n", thresholdValue);
		thresholdValue = newThresholdValue;
		//printf("thresholdValue - After: %d\n", thresholdValue);
		//printf("+++++++++++++++++++\n");

		//calculateThresholdCUDA<<<minMemoryRequired, 32>>>(d_imageArray, thresholdValue, d_greaterThanSum, d_greaterThanCount, d_lessThanSum);
		//printf("minMemoryRequiredForThresholding: %d\n", minMemoryRequiredForThresholding);
		calculateThresholdCUDA<<<minMemoryRequiredForThresholding, THREADS_PER_BLOCK>>>(d_imageArray, thresholdValue, d_greaterThanSum, d_greaterThanCount, d_lessThanSum);
			//printf("calculateThresholdCUDA Return: %s\n", cudaGetErrorString( cudaGetLastError() ));

		//CUSTOM_CUDA_SAFE_CALL(cudaMemcpy(h_greaterThanSum, d_greaterThanSum, (sizeof(unsigned int) * minMemoryRequiredForThresholding), cudaMemcpyDeviceToHost));
		//CUSTOM_CUDA_SAFE_CALL(cudaMemcpy(h_greaterThanCount, d_greaterThanCount, (sizeof(unsigned int) * minMemoryRequiredForThresholding), cudaMemcpyDeviceToHost));
		//CUSTOM_CUDA_SAFE_CALL(cudaMemcpy(h_lessThanSum, d_lessThanSum, (sizeof(unsigned int) * minMemoryRequiredForThresholding), cudaMemcpyDeviceToHost));
		//
		//unsigned int backgroundSetValue, backgroundSetCount, objectSetValue;
		//backgroundSetValue = backgroundSetCount = objectSetValue = 0;
		//for (int i=0 ; i < minMemoryRequiredForThresholding; i++)
		//{
		//	backgroundSetValue += h_greaterThanSum[i];
		//	backgroundSetCount += h_greaterThanCount[i];
		//	objectSetValue += h_lessThanSum[i];
		//}

		//printf("------------RESULTS---------------\n");
		//printf("minMemoryRequired: %d\n", minMemoryRequiredForThresholding);
		//printf("objectSetValue: %d\n", objectSetValue);
		//printf("objectSetCount: %d\n", (307200 - backgroundSetCount));
		//
		//printf("backgroundSetValue: %d\n", backgroundSetValue);
		//printf("backgroundSetCount: %d\n", backgroundSetCount);
		//printf("\n");

		//printf("minBlocksRequiredForThresholdingAtStage2: %d\n", minBlocksRequiredForThresholdingAtStage2);
		calculateThresholdCUDAStage2<<<minBlocksRequiredForThresholdingAtStage2, THREADS_PER_BLOCK>>>(d_greaterThanSum, d_greaterThanCount, d_lessThanSum, \
																									d_greaterThanSumS2, d_greaterThanCountS2, d_lessThanSumS2);


		//CUSTOM_CUDA_SAFE_CALL(cudaMemcpy(h_greaterThanSum, d_greaterThanSumS2, (sizeof(unsigned int) * minBlocksRequiredForThresholdingAtStage2), cudaMemcpyDeviceToHost));
		//CUSTOM_CUDA_SAFE_CALL(cudaMemcpy(h_greaterThanCount, d_greaterThanCountS2, (sizeof(unsigned int) * minBlocksRequiredForThresholdingAtStage2), cudaMemcpyDeviceToHost));
		//CUSTOM_CUDA_SAFE_CALL(cudaMemcpy(h_lessThanSum, d_lessThanSumS2, (sizeof(unsigned int) * minBlocksRequiredForThresholdingAtStage2), cudaMemcpyDeviceToHost));
		//
		//unsigned int backgroundSetValue, backgroundSetCount, objectSetValue;
		//backgroundSetValue = backgroundSetCount = objectSetValue = 0;
		//for (int i=0 ; i < minBlocksRequiredForThresholdingAtStage2; i++)
		//{
		//	backgroundSetValue += h_greaterThanSum[i];
		//	backgroundSetCount += h_greaterThanCount[i];
		//	objectSetValue += h_lessThanSum[i];
		//}

		//printf("------------RESULTS---------------\n");
		//printf("minMemoryRequired: %d\n", minBlocksRequiredForThresholdingAtStage2);
		//printf("objectSetValue: %d\n", objectSetValue);
		//printf("objectSetCount: %d\n", (307200 - backgroundSetCount));
		//
		//printf("backgroundSetValue: %d\n", backgroundSetValue);
		//printf("backgroundSetCount: %d\n", backgroundSetCount);
		//printf("\n");

		//printf("minSharedMemoryRequired * sizeof(int): %d\n", sharedMemSize);
		calculateThresholdCUDAFinal<<<1, 1, sharedMemSize>>>(d_greaterThanSumS2, d_greaterThanCountS2, d_lessThanSumS2, minBlocksRequiredForThresholdingAtStage2, d_threshold);
			//printf("calculateThresholdCUDA Return: %s\n", cudaGetErrorString( cudaGetLastError() ));
		CUSTOM_CUDA_SAFE_CALL(cudaMemcpy(h_threshold, d_threshold, (sizeof(unsigned char)), cudaMemcpyDeviceToHost));

		newThresholdValue = *h_threshold;

		//CUSTOM_CUDA_SAFE_CALL(cudaMemcpy(h_greaterThanSum, d_greaterThanSum, (sizeof(unsigned int) * minMemoryRequiredForThresholding), cudaMemcpyDeviceToHost));
		//CUSTOM_CUDA_SAFE_CALL(cudaMemcpy(h_greaterThanCount, d_greaterThanCount, (sizeof(unsigned int) * minMemoryRequiredForThresholding), cudaMemcpyDeviceToHost));
		//CUSTOM_CUDA_SAFE_CALL(cudaMemcpy(h_lessThanSum, d_lessThanSum, (sizeof(unsigned int) * minMemoryRequiredForThresholding), cudaMemcpyDeviceToHost));


		//for (i = 0; i < MAX_DATA_SIZE; i++)
		//{
		//	if (imageArray[i] > thresholdValue)
		//	{
		//		backgroundSetSize++;
		//		backgroundSetValue += imageArray[i];
		//	}
		//	else
		//	{
		//		objectSetSize++;
		//		objectSetValue += imageArray[i];
		//	}
		//}

		//printf("------------RESULTS---------------\n");
		//printf("minMemoryRequired: %d\n", minMemoryRequired);
		//printf("objectSetValue: %d\n", objectSetValue);
		//printf("objectSetCount: %d\n", (307200 - backgroundSetCount));
		//
		//printf("backgroundSetValue: %d\n", backgroundSetValue);
		//printf("backgroundSetCount: %d\n", backgroundSetCount);
		//printf("\n");

		
		////backgroundSetCount = (backgroundSetCount == 0) ? 1 : backgroundSetCount;
		//backgroundSetAverage = backgroundSetValue / backgroundSetCount;
		////printf("backgroundSetAverage: %d\n", backgroundSetAverage);
		//objectSetAverage = objectSetValue / (MAX_DATA_SIZE - backgroundSetCount);
		////printf("objectSetAverage: %d\n", objectSetAverage);

		//newThresholdValue = (unsigned char) ((objectSetAverage + backgroundSetAverage) / 2);
		//printf("thresholdValue: %d\n", thresholdValue);
		//printf("newThresholdValue: %d\n", newThresholdValue);
		//printf("----------------\n");

		if (thresholdValue == newThresholdValue)
		{
			thresholdFound = 1;
			//printf("\nnewThresholdValue: %d\n",newThresholdValue);
		}

	}

	CUSTOM_CUDA_SAFE_CALL(cudaFree(d_greaterThanSum));
	CUSTOM_CUDA_SAFE_CALL(cudaFree(d_greaterThanCount));
	CUSTOM_CUDA_SAFE_CALL(cudaFree(d_lessThanSum));

	CUSTOM_CUDA_SAFE_CALL(cudaFree(d_greaterThanSumS2));
	CUSTOM_CUDA_SAFE_CALL(cudaFree(d_greaterThanCountS2));
	CUSTOM_CUDA_SAFE_CALL(cudaFree(d_lessThanSumS2));
	
	free(h_threshold);
	CUSTOM_CUDA_SAFE_CALL(cudaFree(d_threshold));

	return (int) newThresholdValue;
}


int iterativeCalculateThresholdUsingArray(unsigned char *imageArray)
{
	int i, thresholdValue, newThresholdValue, objectSetSize, backgroundSetSize, objectSetValue, backgroundSetValue,
		objectSetAverage, backgroundSetAverage, thresholdFound;

	thresholdValue = -1;
	newThresholdValue = 128;  //Half the dynamic range
	thresholdFound = 0;

	//int initialAvgGrayLevel = 0;
	//for (int j = 0; j < (MAX_DATA_SIZE); j++)
	//	initialAvgGrayLevel += imageArray[j];

	//printf("TotalGrayLevel: %d\n", initialAvgGrayLevel);
	//printf("MAX_DATA_SIZE: %d\n", MAX_DATA_SIZE);
	//printf("initialAvgGrayLevel: %d\n", (initialAvgGrayLevel / (MAX_DATA_SIZE)));

	//newThresholdValue = (initialAvgGrayLevel / (MAX_DATA_SIZE));


	while(thresholdFound != 1)
	{
		//printf("\n****************\n");
		//printf("thresholdValue - Before: %d\n", thresholdValue);
		thresholdValue = newThresholdValue;
		//printf("thresholdValue - After: %d\n", thresholdValue);
		
		objectSetSize = 0;
		backgroundSetSize = 0;
		objectSetValue = 0;
		backgroundSetValue = 0;
		
		for (i = 0; i < MAX_DATA_SIZE; i++)
		{
			if (imageArray[i] > thresholdValue)
			{
				backgroundSetSize++;
				backgroundSetValue += imageArray[i];
			}
			else
			{
				objectSetSize++;
				objectSetValue += imageArray[i];
			}
		}

		//printf("objectSetSize: %d\n", objectSetSize);
		//printf("objectSetValue: %d\n", objectSetValue);
		//printf("backgroundSetSize: %d\n", backgroundSetSize);
		//printf("backgroundSetValue: %d\n", backgroundSetValue);
		//printf("\n");

		// Get average of background and object sets
		objectSetAverage = objectSetValue / objectSetSize;
		//printf("objectSetAverage: %d\n", objectSetAverage);

		//backgroundSetSize = (backgroundSetSize == 0) ? 1 : backgroundSetSize;
		backgroundSetAverage = backgroundSetValue / backgroundSetSize;
		//printf("backgroundSetAverage: %d\n", backgroundSetAverage);

		newThresholdValue = (objectSetAverage + backgroundSetAverage) / 2;
		//printf("newThresholdValue: %d\n",newThresholdValue);
		//printf("----------------\n");

		if (thresholdValue == newThresholdValue)
		{
			thresholdFound = 1;
			//printf("\nnewThresholdValue: %d\n",newThresholdValue);
		}

	}

	return newThresholdValue;
}

// MAX_DATA_SIZE = 640*480
// THREADS_PER_BLOCK = 256
// BLOCKS_PER_GRID_ROW (MAX_DATA_SIZE/THREADS_PER_BLOCK) = 1200
// segmentImageCoalescedCUDA<<<BLOCKS_PER_GRID_ROW, THREADS_PER_BLOCK>>>(d_imageArray, threshold);
__global__ void segmentImageCUDA(unsigned char *d_imageArray, unsigned char threshold)
{
	int idx = (blockIdx.x * blockDim.x) + threadIdx.x;

	if (d_imageArray[idx] > threshold)
		d_imageArray[idx] = 0;
	else
		d_imageArray[idx] = 1;
}

__global__ void segmentImageCoalescedCUDA(uchar4 *d_imageArray, unsigned char threshold)
{
	int idx = (blockIdx.x * blockDim.x) + threadIdx.x;

	uchar4 data = d_imageArray[idx];

	data.x = (data.x > threshold ? 0 : 1);
	data.y = (data.y > threshold ? 0 : 1);
	data.z = (data.z > threshold ? 0 : 1);
	data.w = (data.w > threshold ? 0 : 1);

	d_imageArray[idx] = data;
}

__global__ void segmentImageCoalescedCUDATest(unsigned char *d_imageArray, unsigned char threshold)
{
	int idx = (blockIdx.x * blockDim.x) + threadIdx.x;


				//printf("threshold: %d\n", threshold);
				//printf("d_imageArray[0]: %d\n", *(d_imageArray + idx));
				//printf("d_imageArray[1]: %d\n", *(d_imageArray + idx + 1));
				//printf("d_imageArray[2]: %d\n", *(d_imageArray + idx + 2));
				//printf("d_imageArray[3]: %d\n", *(d_imageArray + idx + 3));

				//printf("idx: %d\n", idx);
				//printf("d_imageArray+idx: %d\n", d_imageArray+idx);
				//printf("d_imageArray+(idx*4): %d\n", d_imageArray+(idx*4));

				//unsigned int* iptr = (unsigned int*)(d_imageArray+(idx*4));
				//unsigned int i = *iptr;
				//unsigned char* cptr = (unsigned char*) &i;

	unsigned int* iptr = (unsigned int*)(d_imageArray+(idx*4));
	unsigned int i = *iptr;

	//unsigned char* cptr = (unsigned char*) &i;
	//unsigned char cArray[4];
	typedef __align__(4) unsigned char alignedArray;
	alignedArray cArray[4];

	unsigned int mask = 0x000000FF;

	//cArray[0] = (unsigned char) (i & mask);
	//for(int s=0; s<4; s++)
	//for(unsigned char s=0; s<4; s++)
	//{
	//	if (s==0)
	//		cArray[s] = (unsigned char) (i & mask);
	//	else
	//	{
	//		mask = mask << 8;
	//		cArray[s] = (unsigned char) ((i & mask) >> (8*s));
	//	}
	//}

	//uchar1  sresult;
	//sresult = make_uchar1(0);
	//for(sresult.x=0; sresult.x<4; sresult.x++)
	//{
	//		cArray[sresult.x] = ((i & mask) >> (8*sresult.x));
	//		mask = mask << 8;
	//}


	// Unaligned access
	//for(unsigned char s=0; s<4; s++)
	//{
	//		cArray[s] = (unsigned char) ((i & mask) >> (8*s));
	//		mask = mask << 8;
	//}


	// Unaligned access - however the compiler seems to correct the memory address access by masking and shifting
	// addr[n] replaced with ((*(addr & ~3)) >> (addr & 3)) & 0xff; -- from irc

	//mask = mask << 0;
	//cArray[0] = (unsigned char) (i & mask) >> 0;
	//mask = mask << 8;
	//cArray[1] = (unsigned char) ((i & mask) >> 8);
	//mask = mask << 8;
	//cArray[2] = (unsigned char) ((i & mask) >> 16);
	//mask = mask << 8;
	//cArray[3] = (unsigned char) ((i & mask) >> 24);

	//printf("i: %ud\n", i);

    for(int s=0; s<4; s++) 
    {
		//printf("cArray[%d] before: %d\n", s, cArray[s]);
		if (cArray[s] > threshold)
			cArray[s] = 0;
		else
			cArray[s] = 1;
		//printf("cArray[%d] after : %d\n", s, cArray[s]);
    }

	//printf("i: %ud\n", i);
	//printf("---------------------\n");
	unsigned int* cptr = (unsigned int*) cArray;
	*iptr = *cptr;
	//*(d_imageArray+idx) = *iptr;
}

/*	Cannot coalesce memory reads for unsigned char's as far as I know

	//// Image Size = 640*480
	//// Threads per block = 256
	//// Blocks per grid = 300
	//// segmentImageCoalescedCUDA<<<300, THREADS_PER_BLOCK>>>(d_imageArray, threshold);
	//__global__ void segmentImageCoalescedCUDA(unsigned char *d_imageArray, unsigned char threshold)
	//{
	//	unsigned char temp[4];
	//
	//	int idx = (blockIdx.x * blockDim.x) + threadIdx.x;
	//
	//    for(int s=0; s<4; s++) 
	//    {
	//		if (d_imageArray[(4*idx)+s] > threshold)
	//			temp[s] = 0;
	//		else
	//			temp[s] = 1;
	//    }
	//
	//	//__syncthreads();
	//
	//	for(int s=0; s<4; s++) 
	//		d_imageArray[(4*idx)+s] = temp[s];
	//}
*/



////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* Sequential Methods																																		  */
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


// Create a histogram of the image
short* createHistogram(t_bmp *img)
{
	short *moyLevel = (short*)malloc(sizeof(short) * 256);

	int i, j;

	// Initialise moy array
	for (i = 0; i < 256; i++)
	{
		moyLevel[i] = 0;
	}

	// Get total number of pixels for each moy level
	for (i = 0; i < img->height; i++)
	{
		for (j = 0; j < img->width; j++)
		{
			moyLevel[img->data[i][j].moy]++;
		}
	}

	return moyLevel;
}

// Contrast Enhancement - Linear Stretching
void enhanceContrast(t_bmp *img, short *histogram)
{
	//int outputValue;
	int outputLut;
	int i;
	float startPoint, endPoint, gradient, intercept;
	startPoint = -1;
	endPoint = -1;

	// Calculate I1 - starting point of histogram
	i = 0;
	while(startPoint == -1)
	{
		if (histogram[i] != 0)
		{
			startPoint = (float) i;
		}
		i++;
	}

	// Calculate I2 - end point of histogram
	i = 255;
	while(endPoint == -1)
	{
		if (histogram[i] != 0)
		{
			endPoint = (float) i;
		}
		i--;
	}
	
	//printf("start of histogram: %f \n", startPoint);
	//printf("end of histogram: %f \n", endPoint);

	// M = gradient 
	gradient = 255 / (endPoint - startPoint);
	// C = Y axis intercept
	intercept = -gradient * startPoint;

	//printf("Gradient Standard: %f \n", gradient);
	//printf("Intercept Standard: %f \n", intercept);

	int lut[256];
	for(int i=0; i<256; i++)
		lut[i] = (int) ((gradient * i) + intercept);

	// Modify the image using new values O = gradient * input + intercept
	// Contract enhancement using linear stretching
	for(int i = 0; i < img->height; i++)
	{
		for(int j = 0; j < img->width; j++)
		{
			//outputValue = (int) ((gradient * img->data[i][j].moy) + intercept);
			outputLut = lut[img->data[i][j].moy];
			if (img->data[i][j].moy < (-intercept)/gradient)
			{
				img->data[i][j].moy = 0;
				img->data[i][j].r = 0;
				img->data[i][j].g = 0;
				img->data[i][j].b = 0;
			}
			else if (img->data[i][j].moy > (255 - intercept)/gradient)
			{
				img->data[i][j].moy = 255;
				img->data[i][j].r = 255;
				img->data[i][j].g = 255;
				img->data[i][j].b = 255;
			}
			else
			{
				img->data[i][j].moy = outputLut;
				//// Don't do this - this will generate a bug because we are assigning r to the new value of moy instead of the old one.
				//img->data[i][j].moy = lut[img->data[i][j].moy];
				//img->data[i][j].r = lut[img->data[i][j].moy]
				img->data[i][j].r = outputLut;
				img->data[i][j].g = outputLut;
				img->data[i][j].b = outputLut;

				//img->data[i][j].moy = outputValue;
				//img->data[i][j].r = outputValue;
				//img->data[i][j].g = outputValue;
				//img->data[i][j].b = outputValue;				
			}
		}
	}
}

/// <summary>
/// Carries out a low pass filter on the given image.
/// </summary>
void lowPassFilterImage(t_bmp *img, t_bmp *outputImg)
{
	// Mask size of 9.
	short filter[9];
	

	int count;

	for (int i = 0; i < img->height; i++)
	{
		for (int j = 0; j < img->width; j++)
		{
			count = 0;
			memset(filter, 0, sizeof(short)*9);

			if (i != 479)				 { filter[0] = img->data[i+1][j].moy; count++; }	// Pixel above
			if (i != 0)					 { filter[1] = img->data[i-1][j].moy; count++; }	// Pixel below
			if (j != 0)					 { filter[2] = img->data[i][j-1].moy; count++; }	// Left pixel
			if (j != 639)				 { filter[3] = img->data[i][j+1].moy; count++; }	// Right pixel

			if ((i != 0) && (j != 0))	 { filter[5] = img->data[i-1][j-1].moy; count++; }	// Bottom left
			if ((i != 0) && (j != 639))	 { filter[6] = img->data[i-1][j+1].moy; count++; }	// Bottom pixel		
			if ((i != 479) && (j != 0))	 { filter[7] = img->data[i+1][j-1].moy;	count++; }	// Top Left
			if ((i != 479) && (j != 639)){ filter[8] = img->data[i+1][j+1].moy; count++; }	// Top Right

			filter[4] = img->data[i][j].moy; count++;		// Middle

			short lowPassValue = (filter[0] + filter[1] + filter[2] + filter[3] + filter[4] + filter[5] + filter[6] + filter[7] + filter[8])/count;

			//printf("i: %d - j: %d\n", i, j);
			//printf("Top Left    : %d\n", (filter[7]));
			//printf("Top         : %d\n", (filter[0]));
			//printf("Top Right   : %d\n", (filter[8]));
			//printf("    Left    : %d\n", (filter[2]));

			//printf("    Middle  : %d\n", (filter[4]));
			//printf("    Right   : %d\n", (filter[3]));
			//printf("Bottom Left : %d\n", (filter[5]));
			//printf("Bottom      : %d\n", (filter[1]));
			//printf("Bottom Right: %d\n", (filter[6]));
			//printf("lowPassValue Sum: %d\n", (filter[0] + filter[1] + filter[2] + filter[3] + filter[4] + filter[5] + filter[6] + filter[7] + filter[8]));
			//printf("count       : %d\n", count);
			//printf("lowPassValue: %d\n", lowPassValue);
			//printf("\n");

			//lowPassValue = 0;

			//if((i==0 && j==639) || (i==0 && j==0))
			//{
			//	outputImg->data[i][j].moy = 255;
			//	outputImg->data[i][j].r = 255;
			//	outputImg->data[i][j].g = 255;
			//	outputImg->data[i][j].b = 255;

			//	libbmp_write("imageStage2 - I.bmp", outputImg);
			//}
			//else
			//{

			outputImg->data[i][j].moy = (unsigned char)lowPassValue;
			outputImg->data[i][j].r = (unsigned char)lowPassValue;
			outputImg->data[i][j].g = (unsigned char)lowPassValue;
			outputImg->data[i][j].b = (unsigned char)lowPassValue;

			//libbmp_write("imageStage2 - I.bmp", outputImg);
			//}
		}
	}
}

/// <summary>
/// Carries out a low pass filter on the given image.
/// </summary>
//void lowPassFilterImage(t_bmp *img)  //This implementation is wrong because it writes the result to the same buffer/image instead of a new one.
//{
//	// Mask size of 9.
//	short filter[9];
//
//	for (int i = 1; i < (img->height - 1); i++)
//	{
//		for (int j = 1; j < (img->width - 1); j++)
//		{
//			filter[0] = img->data[i-1][j-1].moy;
//			filter[1] = img->data[i-1][j].moy;
//			filter[2] = img->data[i-1][j+1].moy;
//			//filter[3] = img->data[i][j-1].moy;
//
//			filter[3] = img->data[1][1].moy;
//			filter[4] = img->data[i][j].moy;
//			filter[5] = img->data[i][j+1].moy;
//			filter[6] = img->data[i+1][j-1].moy;
//			filter[7] = img->data[i+1][j].moy;
//			filter[8] = img->data[i+1][j+1].moy;
//
//			short lowPassValue = 0;
//			lowPassValue = (filter[0] + filter[1] + filter[2] + filter[3] + filter[4]
//								  + filter[5] + filter[6] + filter[7] + filter[8])/9;
//
//			if(i==2 && j==638)
//			{
//				img->data[i][j].moy = 0;
//				img->data[i][j].r = 0;
//				img->data[i][j].g = 0;
//				img->data[i][j].b = 255;
//			}
//
//			else if(i==2 && j==639)
//			{
//				printf("Here");
//			}
//
//			else if (i==1 && j==2)
//			{
//				printf("Top Left    : %d\n", (filter[0]));
//				printf("Top         : %d\n", (filter[1]));
//				printf("Top Right   : %d\n", (filter[2]));
//				printf("    Left    : %d\n", (filter[3]));
//				printf("    Middle  : %d\n", (filter[4]));
//				printf("    Right   : %d\n", (filter[5]));
//				printf("Bottom Left : %d\n", (filter[6]));
//				printf("Bottom      : %d\n", (filter[7]));
//				printf("Bottom Right: %d\n", (filter[8]));
//				printf("lowPassValue Sum: %d\n", (filter[0] + filter[1] + filter[2] + filter[3] + filter[4] + filter[5] + filter[6] + filter[7] + filter[8]));
//				printf("lowPassValue: %d\n", lowPassValue);
//				printf("\n");
//
//				img->data[i][j].moy = 0;
//				img->data[i][j].r = 0;
//				img->data[i][j].g = 0;
//				img->data[i][j].b = 0;
//			}
//			else
//			{
//				img->data[i][j].moy = (unsigned char)lowPassValue;
//				img->data[i][j].r = (unsigned char)lowPassValue;
//				img->data[i][j].g = (unsigned char)lowPassValue;
//				img->data[i][j].b = (unsigned char)lowPassValue;
//			}
//		}
//	}
//}

// Median Filter the image based upon moy values
void medianFilterImage(t_bmp *img)
{
	int i, j;
	//int debugX;
	short filter[9];

	for (i = 1; i < (img->height - 1); i++)
	{
		for (j = 1; j < (img->width - 1); j++)
		{ 
			filter[0] = img->data[i-1][j-1].moy;
			filter[1] = img->data[i-1][j].moy;
			filter[2] = img->data[i-1][j+1].moy;
			filter[3] = img->data[i][j-1].moy;
			filter[4] = img->data[i][j].moy;
			filter[5] = img->data[i][j+1].moy;
			filter[6] = img->data[i+1][j-1].moy;
			filter[7] = img->data[i+1][j].moy;
			filter[8] = img->data[i+1][j+1].moy;

			////////debug
			/*
			printf("debug info unsorted \n");
			for(debugX = 0; debugX<9; debugX++)
			{

				printf("value %d is %d \n", debugX, filter[debugX]);
			}*/
			/////////////

			sortArray(filter, 9);

			////////debug
			/*
			printf("debug info sorted\n");
			for(debugX = 0; debugX<9; debugX++)
			{

				printf("value %d is %d \n", debugX, filter[debugX]);
			}*/
			/////////////

			img->data[i][j].moy = (unsigned char)filter[4];
			img->data[i][j].r = (unsigned char)filter[4];
			img->data[i][j].g = (unsigned char)filter[4];
			img->data[i][j].b = (unsigned char)filter[4];

		}
	}
}

// Calculate an automatic threshold
int iterativeCalculateThreshold(t_bmp *img)
{
	int i, j, thresholdValue, newThresholdValue, objectSetSize, backgroundSetSize, objectSetValue,
		backgroundSetValue, objectSetAverage, backgroundSetAverage, thresholdFound;
	newThresholdValue = 100;
	thresholdFound = 0;

	while(thresholdFound != 1)
	{
		thresholdValue = newThresholdValue;
		objectSetSize = 0;
		backgroundSetSize = 0;
		objectSetValue = 0;
		backgroundSetValue = 0;

		for (i = 0; i < img->height; i++)
		{
			// Get size and value of background and object sets
			for (j = 0; j < img->width; j++)
			{
				if (img->data[i][j].moy > thresholdValue)
				{
					objectSetSize++;
					objectSetValue += img->data[i][j].moy;
				}
				else
				{
					backgroundSetSize++;
					backgroundSetValue += img->data[i][j].moy;
				}
			}
		}

		// Get average of background and object sets
		objectSetAverage = objectSetValue / objectSetSize;
		backgroundSetAverage = backgroundSetValue / backgroundSetSize;

		newThresholdValue = (objectSetAverage + backgroundSetAverage) / 2;

		if (thresholdValue == newThresholdValue)
		{
			thresholdFound = 1;
		}

	}

	// fail safe
	if (newThresholdValue > 150)
	{
		newThresholdValue = 140;
	}
	return newThresholdValue;
}

// Segment the image
void segmentImage(t_bmp *img)
{
	int i, j;
	unsigned char threshold = iterativeCalculateThreshold(img);
	//unsigned char threshold = 100;

	for (i = 0; i < img->height; i++)
	{
		for (j = 0; j < img->width; j++)
		{
			// sets pixel to white ? -- black
			if (img->data[i][j].moy > threshold)
			{
				img->data[i][j].r = (unsigned char)0;
				img->data[i][j].g = (unsigned char)0;
				img->data[i][j].b = (unsigned char)0;
				img->data[i][j].moy = (unsigned char)0;
			}
			// sets pixel to black? -- white
			else
			{
				img->data[i][j].r = (unsigned char)255;
				img->data[i][j].g = (unsigned char)255;
				img->data[i][j].b = (unsigned char)255;
				img->data[i][j].moy = (unsigned char)255;
			}
		}
	}
}